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Abstract. Wc consider the problem of finding (possibly noii connected) 
discrete surfaces spanning a finite set of discrete boundary curves in 
the three-dimensional space and minimizing (globally) a discrete energy 
involving mean curvature. Although we consider a fairly general class of 
energies, our main focus is on the Willmore energy, i.e. the total squared 
mean curvature. 

Most works in the literature have been devoted to the approximation of 
a surface evolving by the Willmore flow and, in particular, to the ap- 
proximation of the so-called Willmore surfaces, i.e., the critical points of 
the Willmore energy. Our purpose is to address the delicate task of ap- 
proximating global minimizers of the energy under boundary constraints. 
The main contribution of this work is to translate the nonlinear bound- 
ary value problem into an integer linear program, using a natural for- 
mulation involving pairs of elementary triangles chosen in a pre-specified 
dictionary and allowing self-intersection. 

The reason for such strategy is the well-known existence of algorithms 
that can compute global minimizers of a large class of linear optimizar- 
tion problems, however at a significant computational and memory cost. 
The case of integer linear programming is particularly delicate and usual 
strategies consist in relaxing the integral constraint x € {0, 1} into 
X € [0, 1] which is easier to handle. Our work focuses essentially on 
the connection between the integer linear program and its relaxation. 
We prove that: 

— One cannot guarantee the total unimodularity of the constraint ma- 
trix, which is a sufficient condition for the global solution of the 
relaxed linear program to be always integral, and therefore to be a 
solution of the integer program as well; 

— Furthermore, there are actually experimental evidences that, in some 
cases, solving the relaxed problem yields a fractional solution. 

These facts prove that the problem cannot be tackled with classical linear 
programming solvers, but only with pure integer linear solvers. Neverthe- 
less, due to the very specific structure of the constraint matrix here, we 
strongly believe that it should be possible in the future to design ad-hoc 
integer solvers that yield high-definition approximations to solutions of 
several boundary value problems involving mean curvature, in particular 
the Willmore boundary value problem. 
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1 Introduction 

The Willmore energy of an immersed compact oriented surface f : E ^ M.- 
with boundary dS is defined as 



where H is the mean curvature vector on S, k the geodesic curvature on dS, 
and dA, ds the induced area and length metrics on E, dS. The Willmore en- 
ergy of surfaces with or without boundary plays an important role in geom- 
etry, elastic membranes theory, strings theory, and image processing. Among 
the many concrete optimization problems where the Willmore functional ap- 
pears, let us mention for instance the modeling of biological membranes, the 
design of glasses, and the smoothing of meshed surfaces in computer graphics. 
The Willmore energy is the subject of a long-standing research not only due to 
its relevance to some physical situations but also due to its fundamental prop- 
erty of being conformal invariant, which makes it an interesting substitute to 
the area functional in conformal geometry. Critical points of W with respect to 
interior variations are called Willmore surfaces. They are solutions of the Euler- 
Lagrange equation 5W = whose expression is particularly simple when iV = 3: 
AH + 2H{H^ — K) = 0, being K the Gauss curvature. It is known since Blaschke 
and Thomsen [23] that stereographic projections of compact minimal surfaces in 
S"^ C M"* are always Willmore surfaces in M.^. However, Pinkall exhibited in [22] 
an infinite series of compact embedded Willmore surfaces that are not stereo- 
graphic projections of compact embedded minimal surfaces in S^. Yet Kusner 
conjectured that stereographic projections of Lawson's g-holed tori in S'^ 
should be global minimizers of W among all genus g surfaces. This conjecture is 
still open, except of course for the case g = where the round sphere is known 
to be the unique global minimizer. 

The existence of smooth surfaces that minimize the Willmore energy span- 
ning a given boundary and a conormal field has been proved by Schatzle in |27| . 
Following the notations in |27] , we consider a smooth embedded closed oriented 
curve r C together with a smooth unit normal field np G Np and we 
denote as ztF and ztnr their possible orientations. We assume that there ex- 
ist oriented extensions of ztF, zLnp, that is, there are compact oriented sur- 
faces i7+ C with boundary dS± = ztF and conormal vector field 
co^^ = ztnr on dS±. We also assume that there exists a bounded open set 
B D r such that the set 

{S± oriented extensions of {r,nr), Sj^ connected , 



is not empty. The condition on energy ensures that S+ U is an embedding. 

It follows from ^7\, Corollary 1.2, that the Willmore boundary problem as- 
sociated with (F, nr) in B has a solution, i.e., there exists a compact, oriented. 




S+US- C B, yV{S+ U S-) < 8tt} 
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connected, smooth surface S d B with dS = F, cos — np on dS, and 

W{S) ^ mm{W{S), S smooth, S C B, dS ^ T, co^ = np on dS} 

There have been many contributions to the numerical simulation of Willmore 
surfaces in space dimension N = 3. Among them, Hsu, Kusner and Sullivan 
have tested experimentally in [16| the validity of Kusner's conjecture: starting 
from a triangulated polyhedron in M'^ that is close to a Lawson's surface of 
genus g, they let it evolve by a discrete Willmore flow using Brakke's Surface 
Evolver [6] and check that the solution obtained after convergence is W-stable. 
Recent updates that Brakke brought to its program give now the possibility to 
test the flow with various discrete definitions of the mean curvature. Mayer and 
Simonett [19 introduce a finite difference scheme to approximate axisymmet- 
ric solutions of the Willmore flow. Rusu [26] and Clarenz et al. [8] use a finite 
elements approximation of the flow to compute the evolution of surfaces with 
or without boundary. In both works, position and mean curvature vector are 
taken as independent variables, which is also the case of the contribution by 
Verdera et al. [35], where a triangulated surface with a hole in it is restored 
using the following approach: by the coarea formula, the Willmore energy (actu- 
ally a generalization to other curvature exponents) is replaced with the energy 
of an implicit and smooth representation of the surface, and the mean curvature 
term is replaced by the divergence of an unknown fleld that aims to represent 
the normal field. Droske and Rumpf [3] propose a finite element approach to 
the Willmore fiow but replace the standard fiow equation by its level set for- 
mulation. The contribution of Dziuk jlO| is twofold: it provides a finite element 
approximation to the Willmore fiow with or without boundary conditions that 
can handle as well embedded or immersed surfaces (turning the surface problem 
into a quasi-planar problem), and a consistency result showing the convergence 
of both the discrete surface and the discrete Willmore energy to the continuous 
surface and its energy when the approximated surface has enough regularity. 
Bobenko and Schroder ^ use a difference strategy: they introduce a discrete 
notion of mean curvature for triangulated surfaces computed from the circles 
circumscribed to each triangle that shares with the continuous definition a few 
properties, in particular the invariance with respect to the full Mobius group in 
M^. This discrete definition is vertex-based and a discrete fiow can be derived. 
Based also on several axiomatic constraints but using a finite elements frame- 
work, Wardetzky et al. |33j introduce an edge-based discrete Willmore energy 
for triangulated surfaces. Olischlager and Rumpf [5T] introduce a two step time 
discretization of the Willmore fiow that extends to the Willmore case, at least 
formally, the discrete time approximation of the mean curvature motion due 
to Almgren, Taylor, and Wang [2], and Luckhaus and Sturzenhecker [18]. The 
strategy consists in using the mean curvature fiow to compute an approximation 
of the mean curvature and plug it in a time discrete approximation of the Will- 
more flow. Grzibovskis and Heintz [TJ] , and Esedoglu et al. [TT] discuss how 4th 
order flows can be approximated by iterative convolution with suitable kernels 
and thresholding. 
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While all the previous approaches yield approximations of critical points 
of the Willmore energy, our motivation in this paper is to approximate global 
minimizers of the energy. This is an obviously nontrivial task due to the high 
nonlinearity and nonconvexity of the energy. Yet, for the simpler area functional, 
Sullivan [3T| has shown with a calibration argument that the task of finding min- 
imal surfaces can be turned into a linear problem. Even more, when a discrete 
solution is seeked among surfaces that are union of faces in a cubic grid parti- 
tion of K'^, he proved that the minimization of the linear program is equivalent 
to solving a minimum-cost circulation network flow problem, for which efficient 
codes have been developed by Boykov and Kolmogorov [5, after Ford and Fulk- 
erson |12j . Sullivan [31] did not provide experiments in his paper but this was 
done recently by Grady [13], with applications to the segmentation of medical 
images. 

The linear formulation that we propose here is based on two key ideas: the 
concept of surface continuation constraints that has been pioneered by Sulli- 
van |31j and Grady [13] . and the representation of a triangular surface using 
pairs of triangles. With this representation and a suitable definition of discrete 
mean curvature, we are able to turn into a linear formulation the task of mini- 
mizing discrete representations of any functional of the form 



In the expression of Wip{S), x denotes the space variable, n the normal vector 
field on S and H the mean curvature vector. The linear problem we obtain 
involves integer-valued unknowns and does not seem to admit any simple graph- 
based equivalent. We will therefore discuss whether classical strategies for linear 
optimization can be used. 

The paper is organized as follows: in section[2]we discuss both the chosen rep- 
resentation of surfaces and the definition of discrete mean curvature. In section|3] 
we present a first possible approach yielding a quadratic energy. We present in 
section|4]our linear formulation and discuss whether it can be tackled by classical 
linear optimization techniques. 

2 Discrete framework 

2.1 Triangular meshes from a set of pre-defined triangles 

The equivalence shown by Sullivan between finding minimal surfaces and solving 
a flow problem holds true for discrete surfaces defined as a connected set of cell 
faces in a cellular complex discrete representation of the space. We will consider 
here polyhedral surfaces defined as union of triangles with vertices in (a finite 
subset of) the cubic lattice eZ^ where e = ^ is the resolution scale. Not all 




among discrete immersed surfaces with boundary constraints: 



dS — r, CO — rir on dS. 
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possible triangles are allowed but only those respecting a specified limit on the 
maximal edge length. We assume that each triangle, as well as each triangle edge, 
is represented twice, once for each orientation. We let I denote the collection of 
oriented triangles, N = \I\ its cardinality, and M the number of oriented triangle 
edges. The constrained boundary is given as a contiguous oriented set of triangle 
edges. The orientation of the boundary constrains the spanning surfaces since 
we will allow only spanning triangles whose orientation is compatible. 

In this framework, one can represent a triangular mesh as a binary indicator 
vector X — {0, 1}^ where 1 means that the respective triangle is present in 
the mesh, that it is not. Obviously, not all binary indicator vectors can be 
associated with a triangular surface since the corresponding triangles may not 
be contiguous. However, as discussed by Grady [13 and, in a slightly different 
setting, by Sullivan |30|31j , it is possible to write in a linear form the constraint 
that only binary vectors that correspond to surfaces spanning the given boundary 
are considered. We will see that using the same approach here turns the initial 
boundary value problem into a quadratic program. Another formulation will be 
necessary to get a linear problem. 

2.2 Admissible indicator vectors: a first attempt 



Fig. 1. Incidence of oriented triangles and edges, ei is positively incident to the oriented 
triangle, 62 and 63 are negatively incident, and 64 is not incident to the triangle. 

To define the set of admissible indicator vectors, we first consider a relation- 
ship between oriented triangles and oriented edges which is called incidence: a 
triangle is positive incident to an edge if the edge is one of its borders and the 
two agree in orientation. It is negative incident if the edge is one of its borders, 
but in the opposite orientation. Otherwise it is not incident to the edge. For 
example, the triangle in Figure [T] is positive incident to the edge ei, negative 
incident to 62 and 63 and not incident to 64. 

Being defined as above the set of N oriented triangles and their M oriented 

edges, we introduce the matrix B — .at} whose element bij gives 

je{i',-,M} 

account of the incidence between triangle i and edge j. More precisely 

if edge i is an edge of triangle j with same orientation 
if edge i is an edge of triangle j with opposite orientation 
otherwise 
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The knowledge of which edges are present in the set of prescribed boundary 
segments is expressed as a vector r e {—1, 0, l}*''^ with 

1 if the oriented boundary contains the edge j 

with agreeing orientation 
— 1 if the oriented boundary contains the edge — j 

with opposing orientation 
otherwise 

With these notations set up we can now describe the equation system defining 
that a vector x G {0, 1}^ encodes an oriented triangular mesh with the pre- 
specified oriented boundary. This system has one equation for each edge. If the 
edge is not contained in the given boundary, this equation expresses that, among 
all triangles indicated by x that contain the edge, there are as many triangles 
with same orientation as the edge as triangles with opposite orientation. If the 
edge is contained in the boundary with coherent orientation, there must be 
one more positive incident triangle than negative incident. If it is contained with 
opposite orientation, there is one less positive than negative incident. Altogether 
the constraint for edge j can be expressed as the linear equation 

^ij •^i ~ ''i 

i 

and the entire system as 

Bx^r. (1) 

So far, we did not incorporate the conormal constraint. Actually not all conormal 
constraints are possible, exactly like not all discrete curves can be spanned in 
our framework but only union of edges of dictionary triangles, i.e. the collection 
of triangles defined in the previous section that determine the possible surfaces. 
For the conormal constraint, only the conormal vectors that are tangent to dic- 
tionary triangles sharing an edge with the boundary curve are allowed. Then 
the conormal constraint can be easily plugged into our formulation by simply 
imposing the corresponding triangles to be part of the surface, see Figure [2] and 
by defining accordingly a new boundary indicator vector f. 

Denoting as J' the collection of those additional triangles, the complete con- 
straint reads 

'Bx^r 

We discuss in the next section how discrete mean curvature can be evaluated in 
this framework. 



2.3 Discrete mean curvature on triangular meshes 



The various definitions of discrete mean curvature that have been proposed in the 
literature obviously depend on the chosen discrete representations of surfaces. 
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Fig. 2. The boundary and conormal constraints can be imposed by pre-specifying suit- 
able triangles to be part of the surface. 



Presenting and discussing all possible definitions is out of the scope of the present 
paper. The important thing to know is that there is no fully consistent definition: 
the pointwise convergence of mean curvature cannot be guaranteed in general 
but only in specific situations [15120] . Among the many possible definitions, 
we will use the edge-based one proposed by Polthier ^ for it suits with our 
framework. Recalling that, in the smooth case but also for generalized surfaces 
like varifolds [5^, the first variation of the area can be written in terms of the 
mean curvature, the definition due to Polthier of the mean curvature vector at 
an interior edge e of a simplicial surface reads 

il(e) = |e|cos^A^e (3) 

where |e| is the edge-length, 9e is the dihedral angle between the two triangles 
adjacent to e, and Ne is the angle bisecting unit normal vector, i.e., the unit 
vector collinear to the half sum of the two unit vectors normal to the adjacent 
triangles (see figure [s]). Remark that this formula is a discrete counterpart of 
the continuous H = ki + K2 depending on the principal curvatures, which is 
used in many papers for simplicity as definition of mean curvature. When the 
correct continuous definition H — ^{ki + ^2) is used, the formulas above and 
hereafter should be adapted. The justification of formula ([S]) by Polthier [24125) 




Fig. 3. The edge-based definition of a discrete mean curvature vector due to Polth- 
ier ^24^ depends on the dihedral angle 6^ and the angle bisecting unit normal vector 



is as follows: it is exactly the gradient at any point m G e of the area of the 
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two triangles Ti and T2 adjacent to e, and this gradient does not depend on 
the exact position of to. Indeed, one can subdivide Ti, T2 in four triangles T/, 
i S {1, • • • ,4} having to G e as a vertex and such that Ti = T[ U T2 and 
T2 — T^UT^. The area of each triangle is half the product of the opposite edge's 
length and the height. Therefore, if ei is the positively oriented edge opposite to 
m in the triangle T- and Ji, J2 the rotations in the planes of Ti, T2 by ^, the 
area gradients of T/, i e {1, • • • ,4} at to are | JiCi, 5/162, 5^263, 5/264. The 
sum is the total area gradient of Ti U 12 at to and equals 5(Ji6 + J26), which 
coincides with the formula above. 

As discussed by Wardetsky et al. using the Galerkin theory of approximation, 
this discrete mean curvature is an integrated quantity: it scales as A when each 
space dimension is rescaled by A. A pointwise discrete mean curvature rescaling 
as J is given by (see [33] ) 

FP-(6) = ^COs|iV„ 

where denotes the total area of the two triangles adjacent to e. The factor 
3 comes from the fact that, when the mean curvatures are summed up over all 
edges, the area of each triangle is counted three times, once for each edge. Then 

a discrete counterpart of the energy / f{H) dA is given by 

EAe ,3\e\ Oe j,r . , 

-V.(^cos-iV.). (4) 

edges e 

In particular, the edge-based total squared mean curvature is 

3|ei2 a 



E 



Ae 
edges e 



(cos^)^. (5) 



3 A quadratic program for the minimization of the 
discrete Willmore energy 

Ultimately we are aiming at casting the optimization problem in a form that can 
be handled by standard linear optimization software. Having in mind the frame- 
work described above where a discrete surface spanning the prescribed discrete 
boundary is given as a collection of oriented triangles satisfying equation ([2]) and 
chosen among a pre-specified collection of triangles, a somewhat natural direc- 
tion at first glance seems to be solving a quadratic program. Like in section |2.1[ 
let us indeed denote as [xi) the collection of binary variables associated to the 
"dictionary" of triangles {Ti) and define 

— Cij the common edge to two adjacent triangles Ti and T,-; 

— 9ij the corresponding dihedral angle; 

— Nij the angle bisecting unit normal; 
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— Aij the total area of both triangles. 

Then a continuous energy of the form / (p{x, n, H)dA can be discretized as 

Js 



E 



qtj Xj (6) 



with 



^^V5(ey,iYy, ^^ll^ cos -^Nij) ifi^ j are adjacent 
otherwise 



where tp allows to incorporate dependences on each triangle T^'s position and 
unit normal iVj. In particular, the discrete Willmore energy is 



E 



q^jXi Xj (7) 



with 



ui _ , "ttt^Icos if i 7^ j are adjacent 
lij \ ^Aij 2 

otherwise 

Assuming that the maps tp and (p are positive-valued, both energy matrices 
Q = ilij) and = (qfj) are symmetric matrices in and the mini- 

mization of either Q or Q with boundary constraints turns to be the following 
quadratic program with linear and integrality constraints: 

min {Qx,x) 

X 

such that B X ~ r 

= 1 £ J 
a; £{0,1}^ 

We know of no solution to solve this problem efficiently due to the integrality 
constraint. What is worse, even the relaxed problem where one optimizes over 
X £ [0, 1]^ is very hard to solve: terms of the form XiXj with i ^ j are indefinite, 
so (unless Q has a dominant diagonal) the objective function is a non-convex 
one. 

Moreover, a solution to the relaxed problem would not be of practical use: 
already for the 2D-problem of optimizing curvature energies over curves in the 
plane, the respective quadratic program favors fractional solutions. The relax- 
ation would therefore not be useful for solving the integer program. However, in 
this case Amini et al. [3] showed that one can solve a linear program instead. 
This inspired us for the major contribution of this work: to cast the problem as 
an integer linear program. 
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4 An integer linear programming approach 
4.1 Augmented indicator vectors 

The key idea of the proposed integer hnear program is to consider additional 
indicator vectors. Aside from the indicator variables Xi for basic triangles, one 
now also considers entries Xij corresponding to pairs of adjacent triangles. Such a 
pair is called quadrangle in the following. We will denote x the augmented vector 
{xi, ■ ■ ■ , Xn, ■ ■ ■ , Xij, • • • ) where i ^ j run over all indices of adjacent triangles. 
The cost function can be easily written in a linear form with this augmented 
vector, i.e. it reads 



The major problem to overcome is how to set up a system of constraints that 
guarantees consistency of the augmented vector: the indicator variable Xij for 
the pair of triangles i and j should be 1 if and only if both the variables Xi and 
Xj are 1. Otherwise it should be 0. In addition, one again wants to optimize only 
over indicator vectors that correspond to a triangular mesh. 

To encode this in a linear constraint system, a couple of changes are necessary. 
First of all, we will now have a constraint for each pair of triangle and adjacent 
edge. Secondly, edges are no longer oriented. Still, the set of pre-specified indices 
J implies that the orientation of the border is fixed - we still require that for 
each edge of the boundary an adjacent (oriented) triangle is fixed to constrain 
the conormal information. 

To encode the constraint system we introduce a modified notion of incidence. 
We are no longer interested in incidence of triangles and edges. Instead we now 
consider the incidence of both triangles and quadrangles to pairs of triangles and 
(adjacent) edges. 

For convenience, we define that triangles are positive incident to a pair of 
edge and triangle, whereas all quadrangles are negative incident. 

We propose an incidence matrix where lines correspond to pairs (triangle, 
edge) and columns to either triangles or quadrangles. The entries of this incidence 
matrix are either the incidence of a pair (triangle, edge) with a triangle, defined 
as 



or the incidence of a pair (triangle, edge) with a quadrangle, defined as 




with (see the notations of the previous section) 




rf( (triangle fc,edge e), triangle i) 



1 if i = fc, e is an edge of triangle i 
otherwise 



d( (triangle fc,edge e), quadrangle ij) = 







1 \i i = k OT j = k and i, j 
I otherwise 



share e 
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The columns of this incidence matrix are of two types: either with only O's and 
exactly three 1 (a column corresponding to a triangle T, whose three edges are 
found at lines (T, ei), {T, 62), {T, 63)), or with only O's and exactly two (— l)'s (a 
column corresponding to a quadrangle (Ti, T2) that matches with lines (Ti, 612) 
and (T2,ei2)). 

Again, both the conormal constraints and the boundary edges can be imposed 
by imposing additional triangles indexed by a collection of indices. The general 
constraint has the form 

^d((xfe,e),a;,) + ^ ^((2;^, e), = r^^^^), 

where the right-hand side depends whether the edge e is shared by two triangles 
of the surface (and even several quadrangles in case of self-intersection), or be- 
longs to the new boundary indicated by the additional triangles. If e is an inner 
edge, then the sum must be zero due to our definition of d, otherwise there is 
an adjacent triangle, but no adjacent quadrangle, so the right-hand side should 
be 1: 

, J 1 if fc G e is part of the modified boundary 

(k-e) I Q Qt;herwise 

To sum up, we get the following integer linear program: 

min {w,x) 

X 

such that D X = r' 

xj = 1 e J 
X, €{0,1} Vie {i,...,iV} 

where N is the total number of entries in a;, namely all triangles plus all pairs 
of adjacent triangles. It is worth noticing that such formulation allows triangle 
surfaces with self-intersection. 



(8) 



4.2 On the linear programming relaxation 

Solving integer linear programs is an NP-complete problem, see e.g. [2S1 chapter 
18.1]. This implies that, to the noticeable exception of a few particular prob- 
lems [28], no efficient solutions are known. As a consequence one often resorts 
to solving the corresponding linear programming (LP) relaxation, i.e. one drops 
the integrality constraints. In our case this means to solve the problem: 

min {w,x) (9) 
such that D X ^ r' 

xj ==1 Vj e 

0<x,<l Vi e {1,...,7V} 
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or, equivalently, by suitably augmenting D and r' in order to incorporate the 
second constraint Xj — 1, Vj' G J: 

min {w,x) (10) 

X 

such that Dx = f 

0<i, <1 \/ie{l,...,N} 

There are various algorithms for solving this problem, the most classical being 
the simplex algorithm and several interior point algorithms. Let us now discuss 
the conditions under which these relaxed solutions are also solutions of the orig- 
inal integer linear program. Recalling the basics of LP-relaxation [3S], the set of 
admissible solutions 

P = {x e E^, Dx = f, < X < 1} 

is a polyhedron, i.e. a finite intersection of half-spaces in M^. A classical result 
states that minimizing solutions for the linear objective functions can be seeked 
among the extremal points of P only, i.e. its vertices. Denoting the integral 
envelope of P, that is the convex envelope of P n , another classical result 
states that P has integral vertices only (i.e. vertices with integral coordinates) 
if and only if P = Pg 

Since P = {x e M^, Dx = f,0<x<l}, according to Theorem 19.3 in [25] . 
a sufficient condition for having P = Pe is the property of B being totally 
unimodular, i.e. any square submatrix has determinant either 0, —1 or 1. Under 
this condition, any extremal point of P that is a solution of 

min {w,x) 

bx=f, Xi£[0,l] 

has integral coordinates therefore is a solution of the original integer linear pro- 
gram 

min 

Dx=f, £iG{0,l} 

Theorem 19.3 in [28 mentions an interesting characterization of total unimodu- 
larity due to Paul Camion [7]: a matrix is totally unimodular if, and only if, the 
sum of the entries of every Eulerian square submatrix (i.e. with even rows and 
columns) is divisible by four. 

Unfortunately, we can prove that, as soon as the triangle space is rich enough, 
the incidence matrix D does not satisfy Camion's criterion, therefore is not 
totally unimodular, and neither are the matrices for richer triangles spaces. As 
a consequence, there are choices of the triangle space for which the polyhedron 
P = {x G , Dx = f, < i < 1} may have not only integral vertices, or more 
precisely one cannot guarantee this property thanks to total unimodularity. This 
is summarized in the following theorem. 

Theorem 1. The incidence matrix associated with any triangle space where 
each triangle has a large enough number of adjacent neighbors is not totally 
unimodular. 
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Proof. We show in Figure |4] a configuration and, in Table [T] an associated square 
submatrix of the incidence matrix. The sum of entries over each line and the sum 
over each column are even, though the total sum of the matrix entries is not 
divisible by four. By a result of Camion 7 , the incidence matrix is not totally 
unimodular which yields the conclusion according to [5S][Thm 19.3]. Clearly, 
any triangle space for which this configuration can occur is also associated to an 
incidence matrix that is not totally unimodular. 

It is worth noticing that the previous theorem does not imply that the extremal 
points of the polyhedron P are necessarily not all integral. It only states that 
this cannot be guaranteed as usual by the criterion of total unimodularity. We 
will discuss in the next section what additional informations about integrality 
can be obtained from a few experiments that we have done using classical solvers 
for addressing the relaxed linear problem. 




Fig. 4. A configuration in a triangle space with sufficient resolution. The associated 
incidence matrix is Eulerian (see text) but does not satisfy Camion's criterion, thus is 
not totally unimodular. 
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4.3 Testing the relaxed linear problem 

We have tested the relaxed formulation on a few examples at low-resolution 
using the dual simplex method implemented in the CLP solver. The main reason 
for using low-resolution is that the number of triangles becomes significantly 
important as the resolution increases, and both the computational cost and 
the memory requirements tend to become large. Another reason for working 
at low-resolution is that there is no need to go high before finding a case of 
non-integrality. Indeed, consider the examples in figure [Sj integral solutions are 
obtained when the resolution is very low (i.e. when there is no risk to have 
configurations like in figure [4]). In the last configuration, however, the optimal 
solution of the relaxed problem has fractional entries. This confirms that our 
initial problem cannot be addressed though the classical techniques of relaxation, 
and with usual LP solvers. 

4.4 On integer linear programming 

Our results above indicate that, necessarily, integer linear solvers |28|lj should 
be used. These commonly start with solving the linear programming relaxations, 
then derive further valid inequalities (called cuts) and/or apply a branch-and- 
bound scheme. Due to the small number of fractional values that we have ob- 
served in our experiments, it is quite likely that the derivation of a few cuts only 
would give integral solutions. However, we did not test this so far because of 
the running times of this approach: in cases where we get fractional solutions 
the dual simplex method often needs as long as two weeks and up to 12 GB 
memory! From experience with other linear programming problems we consider 
it likely that the interior point methods implemented in commercial solvers will 
be much faster here (we expect less than a day). At the same time, we expect 
the memory consumption to be considerably much higher, so the method would 
most probably be unusable in practice. 

We strongly believe that a specific integer linear solver should be developed 
rather than using general implementations. It is well known that, for a few 
problems like the knapsack problem ^28] [chapter 24.6], their specific structure 
gives rise to ad-hoc efficient approaches. Recalling that our incidence matrix is 
very sparse and well structured (the nonzero entries of each column are either 
exactly two (—1), or exactly three 1) we strongly believe that an efficient integer 
solver can be developed and our approach can be amenable to higher-resolution 
results in the near future. 

5 Conclusion 

We have shown that the minimization under boundary constraints of mean cur- 
vature based energies over surfaces, and in particular the Willmore energy, can 
be cast as an integer linear program. Unfortunately, this integer program is 
not equivalent to its relaxation so the classical LP algorithms offer iro warranty 
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Fig. 5. A series of experiments (the result and the mesh edges) with increasing resolu- 
tion of the triangle space (and various boundary constraints). An integral solution of 
the relaxed problem is obtained by a standard LP-solver in both top cases. As for the 
last case, the triangle space resolution is now large enough for having configurations 
similar to the counterexample of figure [4] And indeed, an optimal solution is found for 
the relaxed problem that is not integral. The mesh on the bottom-right shows actually 
two nested semi-spheres whose triangles have, at least for a few of them, non binary 
labels. 
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that the integer optimal solution will be found. This implies that pure integer 
linear algorithms must be used, which are in general much more involved. We 
believe however that the particular structure of the problem paves the way to a 
dedicated algorithm that would provide high-resolution global minimizers of the 
Willmore boundary problem and generalizations. This is the purpose of future 
research. 
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